Fit VBGE models to back-calculated length-at-age

Load libraries

pkgs <- c("tidyverse", "tidylog", "broom", "RColorBrewer", "ggridges", "viridis",
          "minpack.lm", "patchwork", "ggtext", "brms", "tidybayes", "modelr")

# minpack.lm needed if using nlsLM()
# if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){
# 
#     install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
#   
#   }

invisible(lapply(pkgs, library,character.only = T))

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'tidylog'


The following objects are masked from 'package:dplyr':

    add_count, add_tally, anti_join, count, distinct, distinct_all,
    distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
    full_join, group_by, group_by_all, group_by_at, group_by_if,
    inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
    relocate, rename, rename_all, rename_at, rename_if, rename_with,
    right_join, sample_frac, sample_n, select, select_all, select_at,
    select_if, semi_join, slice, slice_head, slice_max, slice_min,
    slice_sample, slice_tail, summarise, summarise_all, summarise_at,
    summarise_if, summarize, summarize_all, summarize_at, summarize_if,
    tally, top_frac, top_n, transmute, transmute_all, transmute_at,
    transmute_if, ungroup


The following objects are masked from 'package:tidyr':

    drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
    spread, uncount


The following object is masked from 'package:stats':

    filter


Loading required package: viridisLite

Loading required package: Rcpp

Loading 'brms' package (version 2.19.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').


Attaching package: 'brms'


The following object is masked from 'package:stats':

    ar



Attaching package: 'tidybayes'


The following objects are masked from 'package:brms':

    dstudent_t, pstudent_t, qstudent_t, rstudent_t


The following objects are masked from 'package:ggridges':

    scale_point_color_continuous, scale_point_color_discrete,
    scale_point_colour_continuous, scale_point_colour_discrete,
    scale_point_fill_continuous, scale_point_fill_discrete,
    scale_point_size_continuous



Attaching package: 'modelr'


The following object is masked from 'package:broom':

    bootstrap

# devtools::install_github("seananderson/ggsidekick") ## not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Load functions
home <- here::here()
fxn <- list.files(paste0(home, "/R/functions"))
invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))

Load cache

# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/analyze-data/02-fit-vbge_cache/html"))

Read and trim data

d <- #read_csv(paste0(home, "/data/for-analysis/dat.csv")) %>% 
  read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv") %>% 
  filter(age_ring == "Y", # use only length-at-age by filtering on age_ring
         !area %in% c("VN", "TH")) 

Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 43,886 rows (13%), 294,574 rows remaining

# Minimum number of observations per area and cohort
d$area_cohort <- as.factor(paste(d$area, d$cohort))

d <- d %>%
  group_by(area_cohort) %>% 
  filter(n() > 100)

group_by: one grouping variable (area_cohort)
filter (grouped): removed 2,637 rows (1%), 291,937 rows remaining

# Minimum number of observations per area, cohort, and age
d$area_cohort_age <- as.factor(paste(d$area, d$cohort, d$age_bc))

d <- d %>%
  group_by(area_cohort_age) %>% 
  filter(n() > 10)

group_by: one grouping variable (area_cohort_age)
filter (grouped): removed 3,521 rows (1%), 288,416 rows remaining

# Minimum number of cohorts in a given area
cnt <- d %>%
  group_by(area) %>%
  summarise(n = n_distinct(cohort)) %>%
  filter(n >= 10)

group_by: one grouping variable (area)
summarise: now 10 rows and 2 columns, ungrouped
filter: no rows removed

d <- d[d$area %in% cnt$area, ]

# Plot cleaned data
ggplot(d, aes(age_bc, length_mm)) +
  geom_jitter(size = 0.1, height = 0, alpha = 0.1) +
  scale_x_continuous(breaks = seq(20)) +
  theme(axis.text.x = element_text(angle = 0)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  labs(x = "Age", y = "Length (mm)") +
  guides(color = "none") + 
  facet_wrap(~area, scale = "free_x")

# Longitude and latitude attributes for each area
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) %>%
  mutate_at(c("lat", "lon"), as.numeric)

mutate_at: converted 'lat' from character to double (0 new NA)
           converted 'lon' from character to double (0 new NA)

Fit VBGE models

# Get individual growth parameters (functions: VBGF/Gompertz,nls_out,fit_nls)
IVBG <- d %>%
  group_by(ID) %>%
  summarize(nls_out(fit_nls(length_mm, age_bc, min_nage = 5, model = "VBGF")))

group_by: one grouping variable (ID)

summarize: now 75,245 rows and 5 columns, ungrouped

IVBG <- IVBG %>% drop_na(k) # The NAs are because the model didn't converge or because they were below the threshold age

drop_na: removed 51,640 rows (69%), 23,605 rows remaining

IVBG <- IVBG %>%
  mutate(k_lwr = k - 1.96*k_se,
         k_upr = k + 1.96*k_se,
         linf_lwr = linf - 1.96*linf_se,
         linf_upr = linf + 1.96*linf_se,
         row_id = row_number())

mutate: new variable 'k_lwr' (double) with 23,603 unique values and 0% NA

        new variable 'k_upr' (double) with 23,603 unique values and 0% NA

        new variable 'linf_lwr' (double) with 23,603 unique values and 0% NA

        new variable 'linf_upr' (double) with 23,603 unique values and 0% NA

        new variable 'row_id' (integer) with 23,605 unique values and 0% NA

# Plot all K's
IVBG %>%
  #filter(row_id < 5000) %>%
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr)) +
  geom_point(alpha = 0.2) +
  geom_errorbar(alpha = 0.2) +
  NULL

# Plot all L_inf's
IVBG %>%
  #filter(row_id < 5000) %>%
  ggplot(aes(row_id, linf, ymin = linf_lwr, ymax = linf_upr)) +
  geom_point(alpha = 0.2) +
  geom_errorbar(alpha = 0.2) +
  NULL

# We can also consider removing individuals where the SE of k is larger than the fit
IVBG %>%
  #mutate(keep = ifelse(k > quantile(k_se, probs = 0.95), "N", "Y")) %>%
  mutate(keep = ifelse(k < k_se, "N", "Y")) %>%
  #filter(row_id < 10000) %>%
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~keep, ncol = 1) +
  geom_errorbar(alpha = 0.5) +
  NULL

mutate: new variable 'keep' (character) with 2 unique values and 0% NA

# Add back cohort and area variables
IVBG <- IVBG %>% 
  left_join(d %>% select(ID, area_cohort)) %>% 
  separate(area_cohort, into = c("area", "cohort"), remove = TRUE, sep = " ") %>% 
  mutate(cohort = as.numeric(cohort))

Adding missing grouping variables: `area_cohort_age`

select: dropped 11 variables (length_mm, age_bc, age_catch, age_ring, area, …)

Joining with `by = join_by(ID)`
left_join: added 2 columns (area_cohort_age, area_cohort)
> rows only in x 0
> rows only in y (146,393)
> matched rows 142,023 (includes duplicates)
> =========
> rows total 142,023
mutate: converted 'cohort' from character to double (0 new NA)

# Compare how the means and quantiles differ depending on this filtering
IVBG_filter <- IVBG %>%
  drop_na(k_se) %>%
  #filter(k_se < quantile(k_se, probs = 0.95)) %>% 
  filter(k_se < k)

drop_na: no rows removed
filter: removed 4,192 rows (3%), 137,831 rows remaining

# Summarize growth coefficients by cohort and area
VBG <- IVBG %>%
  filter(k_se < k) %>% # new!
  group_by(cohort, area) %>%
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            linf_median = quantile(linf, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) %>% 
  ungroup()

filter: removed 4,192 rows (3%), 137,831 rows remaining
group_by: 2 grouping variables (cohort, area)
summarize: now 306 rows and 7 columns, one group variable remaining (cohort)
ungroup: no grouping variables

VBG_filter <- IVBG_filter %>%
  group_by(cohort, area) %>%
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) %>% 
  ungroup()

group_by: 2 grouping variables (cohort, area)
summarize: now 306 rows and 6 columns, one group variable remaining (cohort)
ungroup: no grouping variables

ggplot() +
  geom_ribbon(data = VBG, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,
                              fill = "All k's"), alpha = 0.4) +
  geom_ribbon(data = VBG_filter, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,
                                     fill = "Filtered"), alpha = 0.4) +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered"), linetype = 2) + 
  guides(fill = "none") +
  facet_wrap(~area)

ggplot() +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered"), linetype = 2) + 
  guides(fill = "none") +
  facet_wrap(~area)

# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.

Calculate sample size

sample_size <- IVBG %>%
  group_by(area) %>% 
  summarise(n_cohort = length(unique(cohort)),
            min_cohort = min(cohort),
            max_cohort = max(cohort),
            n_individuals = length(unique(ID)),
            n_data_points = n())

group_by: one grouping variable (area)

summarise: now 10 rows and 6 columns, ungrouped

sample_size

# A tibble: 10 × 6
   area  n_cohort min_cohort max_cohort n_individuals n_data_points
   <chr>    <int>      <dbl>      <dbl>         <int>         <int>
 1 BS          17       1985       2001          1334          7688
 2 BT          22       1972       1995           961          5371
 3 FB          47       1969       2015          6040         37381
 4 FM          37       1963       1999          5057         32642
 5 HO          29       1982       2015          1150          6210
 6 JM          60       1953       2014          4868         28749
 7 MU          18       1981       1999          1104          6666
 8 RA          14       1985       2003           572          3122
 9 SI_EK       24       1968       2011           658          3571
10 SI_HA       38       1967       2013          1861         10623

sample_size %>%
  ungroup() %>%
  summarise(sum_ind = sum(n_individuals), sum_n = sum(n_data_points))

ungroup: no grouping variables

summarise: now one row and 2 columns, ungrouped

# A tibble: 1 × 2
  sum_ind  sum_n
    <int>  <int>
1   23605 142023

write_csv(sample_size, paste0(home, "/output/sample_size.csv"))

Add GAM-predicted temperature to growth models

pred_temp <- read_csv(paste0(home, "/output/gam_predicted_temps.csv")) %>% 
  rename(cohort = year)

Rows: 380 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): area, model
dbl (2): year, temp

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (cohort)

VBG <- VBG %>%
  left_join(pred_temp, by = c("area", "cohort"))

left_join: added 2 columns (temp, model)
           > rows only in x     0
           > rows only in y  ( 74)
           > matched rows     306
           >                 =====
           > rows total       306

# Save data for map-plot
cohort_sample_size <- IVBG %>%
  group_by(area, cohort) %>% 
  summarise(n = n()) # individuals, not samples!

group_by: 2 grouping variables (area, cohort)
summarise: now 306 rows and 3 columns, one group variable remaining (area)

VBG <- left_join(VBG, cohort_sample_size, by = c("cohort", "area"))

left_join: added one column (n)
           > rows only in x     0
           > rows only in y  (  0)
           > matched rows     306
           >                 =====
           > rows total       306

write_csv(VBG, paste0(home, "/output/vbg.csv"))

# Calculate the plotting order (also used for map plot)
order <- VBG %>%
  ungroup() %>%
  group_by(area) %>%
  summarise(min_temp = min(temp)) %>%
  arrange(desc(min_temp))

ungroup: no grouping variables
group_by: one grouping variable (area)
summarise: now 10 rows and 2 columns, ungrouped

order

# A tibble: 10 × 2
   area  min_temp
   <chr>    <dbl>
 1 SI_HA     7.83
 2 BS        6.22
 3 BT        5.87
 4 FM        5.87
 5 SI_EK     5.49
 6 MU        5.16
 7 FB        5.04
 8 HO        3.99
 9 JM        3.37
10 RA        3.06

write_csv(order, paste0(home, "/output/ranked_temps.csv"))

nareas <- length(unique(order$area)) + 2 # to skip the brightest colors that are hard to read
colors <- colorRampPalette(brewer.pal(name = "RdYlBu", n = 10))(nareas)[-c(6,7)]

Plot VBGE fits

# Sample 30 IDs per area and plot their data and VBGE fits
set.seed(4)
ids <- IVBG %>% distinct(ID, .keep_all = TRUE) %>% group_by(area) %>% slice_sample(n = 30)

distinct: removed 118,418 rows (83%), 23,605 rows remaining

group_by: one grouping variable (area)

slice_sample (grouped): removed 23,305 rows (99%), 300 rows remaining

IVBG2 <- IVBG %>%
  filter(ID %in% ids$ID) %>% 
  distinct(ID, .keep_all = TRUE) %>% 
  select(ID, k, linf)

filter: removed 140,319 rows (99%), 1,704 rows remaining

distinct: removed 1,404 rows (82%), 300 rows remaining

select: dropped 10 variables (k_se, linf_se, k_lwr, k_upr, linf_lwr, …)

d2 <- d %>%
  ungroup() %>%
  filter(ID %in% ids$ID) %>%
  left_join(IVBG2, by = "ID") %>%
  mutate(length_mm_pred = linf*(1-exp(-k*age_bc)))

ungroup: no grouping variables

filter: removed 286,712 rows (99%), 1,704 rows remaining

left_join: added 2 columns (k, linf)

           > rows only in x       0

           > rows only in y  (    0)

           > matched rows     1,704

           >                 =======

           > rows total       1,704

mutate: new variable 'length_mm_pred' (double) with 1,697 unique values and 0% NA

fits_ind <- ggplot(d2, aes(age_bc, length_mm, group = ID, color = ID)) +
  geom_jitter(width = 0.3, height = 0, alpha = 0.5, size = 0.4) +
  geom_line(aes(age_bc, length_mm_pred, group = ID, color = ID),
            inherit.aes = FALSE, alpha = 0.8, linewidth = 0.3) + 
  guides(color = "none") +
  scale_color_viridis(discrete = TRUE, name = "Site", option = "cividis") +
  labs(x = "Age (years)", y = "Length (mm)") +
  facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area), ncol = 5) + 
  NULL

k <- IVBG %>% 
  ggplot(aes(factor(area, order$area), k, fill = factor(area, order$area))) + 
  coord_cartesian(ylim = c(0, 0.7)) +
  geom_violin(alpha = 0.8, color = NA) +  
  geom_boxplot(outlier.color = NA, width = 0.2, alpha = 0.9, fill = NA, size = 0.4) +
  scale_fill_manual(values = colors, name = "Site") +
  guides(fill = "none", color = "none") +
  labs(x = "Site", y = expression(italic(k))) + 
  coord_flip()

Coordinate system already present. Adding new coordinate system, which will
replace the existing one.

linf <- IVBG %>% 
  filter(linf < 2300) %>% 
  filter(linf > 130) %>% 
  ggplot(aes(factor(area, order$area), linf, fill = factor(area, order$area))) + 
  geom_violin(alpha = 0.8, color = NA) +  
  geom_boxplot(outlier.color = NA, width = 0.1, alpha = 0.9, fill = NA, size = 0.4) +
  coord_cartesian(ylim = c(0, 2000)) +
  scale_fill_manual(values = colors, name = "Site") +
  guides(fill = "none", color = "none") +
  labs(x = "", y = expression(paste(italic(L[infinity]), " [mm]"))) +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) + 
  coord_flip()

filter: removed 1,218 rows (1%), 140,805 rows remaining

filter: no rows removed

Coordinate system already present. Adding new coordinate system, which will
replace the existing one.

#fits_ind / (k + linf) + plot_annotation(tag_levels = "A") + plot_layout(heights = c(1, 1.8))

k2 <- k + coord_cartesian(ylim = c(0, 0.65))

Coordinate system already present. Adding new coordinate system, which will
replace the existing one.

fits_ind / k2 + plot_annotation(tag_levels = "A") #+ plot_layout(heights = c(1, 1.8))

ggsave(paste0(home, "/figures/vb_pars_fits.pdf" ), width = 17, height = 18, units = "cm")

Fit Sharpe-Schoolfield model to k

dat <- VBG %>%
  select(k_median, temp, area) %>%
  rename(rate = k_median)

select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)

rename: renamed one variable (rate)

Can we fit a single Sharpe Schoolfield with area-specific parameters with brms?

# Again, here are the data we are fitting:
ggplot(dat, aes(temp, rate, color = factor(area, order$area))) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6)  

hist(dat$rate)

# Here's the equation (from the r package rTPC):
# > sharpeschoolhigh_1981
# function (temp, r_tref, e, eh, th, tref) 
# {
#     tref <- 273.15 + tref
#     k <- 8.62e-05
#     boltzmann.term <- r_tref * exp(e/k * (1/tref - 1/(temp + 
#         273.15)))
#     inactivation.term <- 1/(1 + exp(eh/k * (1/(th + 273.15) - 
#         1/(temp + 273.15))))
#     return(boltzmann.term * inactivation.term)
# }

# Add in fixed parameters
dat$bk <- 8.62e-05
dat$tref <- 8 + 273.15

# (better visualization including bounds further down)
n = 10000
hist(rnorm(0.8, 1, n = n), main = "e", xlim = c(0, 5))

hist(rnorm(2, 1, n = n), main = "eh", xlim = c(0, 7))

hist(rnorm(0.25, 0.5, n = n), main = "r_tref", xlim = c(0, 2.5))

hist(rnorm(12, 5, n = n), main = "th", xlim = c(5, 30))

# These work nicely with the pp_check
prior <- c(prior(normal(0.8, 1), nlpar = "e", lb = 0),
           prior(normal(2, 1), nlpar = "eh", lb = 0),
           prior(normal(0.3, 0.5), nlpar = "rtref", lb = 0),
           prior(normal(10, 5), nlpar = "th", lb = 0))

# Now to a prior predictive check
fit_prior <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1, 
     nl = TRUE),
  data = dat,
  cores = 2,
  chains = 2,
  iter = 1500,
  sample_prior = "only",
  seed = 9,
  prior = prior
)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/Rcpp/1.0.10/e749cae40fa9ef469b6050959517453c/Rcpp/include/"  -I"/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/"  -I"/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/unsupported"  -I"/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/BH/include" -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/src/"  -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/"  -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/RcppParallel/5.1.7/a45594a00f5dbb073d5ec9f48592a08a/RcppParallel/include/"  -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/rstan/2.32.3/4ac5b7639d28cd4fab19baaf46f33c6a/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/Dense:1:
In file included from /Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/Core:88:
/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/Dense:1:
/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

# Global prior predictive check in relation to data. Doesn't loo too informative...
dat %>%
  data_grid(temp = seq_range(temp, n = 51)) %>% 
  ungroup() %>% 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) %>% 
  add_epred_draws(fit_prior) %>% 
  ggplot(aes(temp, y = .epred)) +
  stat_lineribbon(.width = c(0.75), alpha = 0.3, color = "black", fill = "black") + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Site") +
  labs(y = "Expectation of the posterior predictive distribution", x = "Temperature (°C)") +
  NULL

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ggsave(paste0(home, "/figures/supp/prior_predictive_check.pdf" ), width = 15, height = 11, units = "cm")

# Now make sure it converges with real data but no random effects
fit_global <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1,
     nl = TRUE),
  data = dat,
  cores = 4,
  chains = 4,
  iter = 4000, 
  seed = 9,
  sample_prior = "yes",
  family = student,
  prior = prior
  )

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

plot(fit_global)

pp_check(fit_global)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

plot(conditional_effects(fit_global, effect = "temp"), points = TRUE)

# Plot prior vs posterior
post <- fit_global %>%
  as_draws_df() %>%
  dplyr::select(b_rtref_Intercept, b_e_Intercept, b_eh_Intercept, b_th_Intercept) %>% 
  rename(r_tref = b_rtref_Intercept, 
         e = b_e_Intercept,
         eh = b_eh_Intercept,
         th = b_th_Intercept) %>% 
  pivot_longer(everything(), names_to = "parameter") %>% 
  mutate(type = "Posterior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

prior_draws <- fit_global %>%
  as_draws_df() %>%
  dplyr::select(prior_b_rtref, prior_b_e, prior_b_eh, prior_b_th) %>% 
  rename(r_tref = prior_b_rtref, 
         e = prior_b_e,
         eh = prior_b_eh,
         th = prior_b_th) %>% 
  pivot_longer(everything(), names_to = "parameter") %>% 
  mutate(type = "Prior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

dist <- bind_rows(prior_draws, post)

ggplot(dist, aes(value, fill = type)) +
  geom_density(color = NA, alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free") +
  theme(legend.position = c(0.35, 0.9)) +
  labs(fill = "")

Ok, seems to work with priors and everything. They are roughly as broad as can be and still have a fit. Now fit the full model, with random area effects, more iterations and chains!

# Now look at random area effects!
nrow(dat)

[1] 306

# Student-t model with random effects
fitbs <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1 + (1|area),
     nl = TRUE),
  data = dat,
  family = student,
  cores = 4,
  chains = 4,
  iter = 5000,
  sample_prior = "yes",
  save_pars = save_pars(all = TRUE),
  seed = 9,
  prior = prior,
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

print(fitbs, digits = 4)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 306) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~area (Number of levels: 10) 
                    Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS
sd(rtref_Intercept)   0.0447    0.0384   0.0020   0.1434 1.0027     1649
sd(e_Intercept)       0.1823    0.1285   0.0089   0.4883 1.0009     2998
sd(eh_Intercept)      0.3051    0.3025   0.0083   1.0850 1.0031     1911
sd(th_Intercept)      0.8736    0.7670   0.0371   2.8345 1.0023     2302
                    Tail_ESS
sd(rtref_Intercept)     2800
sd(e_Intercept)         3025
sd(eh_Intercept)        4173
sd(th_Intercept)        4057

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
rtref_Intercept   0.4328    0.2079   0.2148   1.0089 1.0023     1418     3286
e_Intercept       0.8500    0.4512   0.1991   1.8886 1.0026     1442     3741
eh_Intercept      1.6155    0.4614   0.6903   2.5314 1.0004     3357     3018
th_Intercept      9.1672    3.8499   2.9851  17.2414 1.0021     1413     3481

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
sigma   0.0397    0.0020   0.0357   0.0435 1.0002    10449     7698
nu     24.6696   13.2991   8.2356  59.1964 1.0003    10512     7702

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Inspect the model

# Check fit

library(performance)

Attaching package: 'performance'

The following objects are masked from 'package:modelr':

    mae, mse, rmse

fitbs

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 306) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~area (Number of levels: 10) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(rtref_Intercept)     0.04      0.04     0.00     0.14 1.00     1649     2800
sd(e_Intercept)         0.18      0.13     0.01     0.49 1.00     2998     3025
sd(eh_Intercept)        0.31      0.30     0.01     1.08 1.00     1911     4173
sd(th_Intercept)        0.87      0.77     0.04     2.83 1.00     2302     4057

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept     0.43      0.21     0.21     1.01 1.00     1418     3286
e_Intercept         0.85      0.45     0.20     1.89 1.00     1442     3741
eh_Intercept        1.62      0.46     0.69     2.53 1.00     3357     3018
th_Intercept        9.17      3.85     2.99    17.24 1.00     1413     3481

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.04      0.00     0.04     0.04 1.00    10449     7698
nu       24.67     13.30     8.24    59.20 1.00    10512     7702

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

r2_bayes(fitbs)

# Bayesian R2 with Compatibility Interval

  Conditional R2: 0.373 (95% CI [0.301, 0.438])

r2_bayes(fit_global)

# Bayesian R2 with Compatibility Interval

  Conditional R2: 0.120 (95% CI [0.057, 0.187])

# qq-plots
p1 <- dat %>%
  add_predicted_draws(fitbs) %>%
  summarise(
    p_residual = mean(.prediction < rate),
    z_residual = qnorm(p_residual),
    .groups = "drop_last"
  ) %>%
  ggplot(aes(sample = z_residual)) +
  geom_qq() +
  geom_abline() + 
  labs(x = "Theoretical", y = "Sample") + 
  ggtitle("Mixed model") +
  theme(aspect.ratio = 1)

summarise: now 306 rows and 8 columns, 5 group variables remaining (rate, temp, area, bk, tref)

p2 <- dat %>%
  add_predicted_draws(fit_global) %>%
  summarise(
    p_residual = mean(.prediction < rate),
    z_residual = qnorm(p_residual),
    .groups = "drop_last"
  ) %>%
  ggplot(aes(sample = z_residual)) +
  geom_qq() +
  geom_abline() + 
  labs(x = "Theoretical", y = "Sample") + 
  ggtitle("Global model") +
  theme(aspect.ratio = 1)

summarise: now 306 rows and 8 columns, 5 group variables remaining (rate, temp, area, bk, tref)

p3 <- pp_check(fitbs) + theme(legend.position = c(0.1, 0.9)) + ggtitle("Mixed model")

Using 10 posterior draws for ppc type 'dens_overlay' by default.

p4 <- pp_check(fit_global) + theme(legend.position = c(0.1, 0.9)) + ggtitle("Global model")

Using 10 posterior draws for ppc type 'dens_overlay' by default.

(p1 + p2) / (p3 + p4)

ggsave(paste0(home, "/figures/supp/sharpe_brms_ppcheck_qq.pdf" ), width = 17, height = 17, units = "cm")

# Plot prior vs posterior
post <- fitbs %>%
  as_draws_df() %>%
  dplyr::select(b_rtref_Intercept, b_e_Intercept, b_eh_Intercept, b_th_Intercept) %>% 
  rename(r_tref = b_rtref_Intercept, 
         e = b_e_Intercept,
         eh = b_eh_Intercept,
         th = b_th_Intercept) %>% 
  pivot_longer(everything(), names_to = "parameter") %>% 
  mutate(type = "Posterior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 10000x4, now 40000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

prior_draws <- fitbs %>%
  as_draws_df() %>%
  dplyr::select(prior_b_rtref, prior_b_e, prior_b_eh, prior_b_th) %>% 
    rename(r_tref = prior_b_rtref, 
         e = prior_b_e,
         eh = prior_b_eh,
         th = prior_b_th) %>% 
  pivot_longer(everything(), names_to = "parameter") %>% 
  mutate(type = "Prior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 10000x4, now 40000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

dist <- bind_rows(prior_draws, post)

ggplot(dist, aes(value, fill = type)) +
  geom_density(color = NA, alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free") + 
  theme(legend.position = c(0.35, 0.9)) +
  labs(fill = "", x = "Value", y = "Density")

ggsave(paste0(home, "/figures/supp/sharpe_brms_prior.pdf" ), width = 17, height = 17, units = "cm")

Sensitivity with respect to priors

# Set a bunch of priors with different sigmas
  # prior1 <- c(prior(normal(0.8, 1), nlpar = "e", lb = 0),
  #             prior(normal(2, 1), nlpar = "eh", lb = 0),
  #             prior(normal(0.3, 0.5), nlpar = "rtref", lb = 0),
  #             prior(normal(10, 5), nlpar = "th", lb = 0))

# For each object, we increase sd by a factor 1.5 from previous
  prior2 <- c(prior(normal(0.8, 1.5), nlpar = "e", lb = 0),
              prior(normal(2, 1.5), nlpar = "eh", lb = 0),
              prior(normal(0.3, 0.75), nlpar = "rtref", lb = 0),
              prior(normal(10, 7.5), nlpar = "th", lb = 0))
  
  prior3 <- c(prior(normal(0.8, 3), nlpar = "e", lb = 0),
              prior(normal(2, 3), nlpar = "eh", lb = 0),
              prior(normal(0.3, 1.125), nlpar = "rtref", lb = 0),
              prior(normal(10, 11.25), nlpar = "th", lb = 0))
  
  prior4 <- c(prior(normal(0.8, 4.5), nlpar = "e", lb = 0),
              prior(normal(2, 4.5), nlpar = "eh", lb = 0),
              prior(normal(0.3, 1.6875), nlpar = "rtref", lb = 0),
              prior(normal(10, 16.875), nlpar = "th", lb = 0))
  
  prior5 <- c(prior(normal(0.8, 6.75), nlpar = "e", lb = 0),
              prior(normal(2, 6.75), nlpar = "eh", lb = 0),
              prior(normal(0.3, 2.53125), nlpar = "rtref", lb = 0),
              prior(normal(10, 25.3125), nlpar = "th", lb = 0))
  
  prior6 <- c(prior(normal(0.8, 10.125), nlpar = "e", lb = 0),
              prior(normal(2, 10.125), nlpar = "eh", lb = 0),
              prior(normal(0.3, 3.796875), nlpar = "rtref", lb = 0),
              prior(normal(10, 37.96875), nlpar = "th", lb = 0))

  priors <- list(prior2, prior3, prior4, prior5, prior6)

# List for storing posteriors and priors
dist_list <- list()
  
for(i in 1:5) {
  
  prior <- priors[[i]]
  
  fit_sens <- brm(
    bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
       rtref + e + eh + th ~ 1,
       nl = TRUE),
    data = dat,
    cores = 4,
    chains = 4,
    iter = 4000, 
    seed = 9,
    sample_prior = "yes",
    family = student,
    prior = prior
    )
  
  # Plot prior vs posterior
  post <- fit_sens %>%
    as_draws_df() %>%
    dplyr::select(b_rtref_Intercept, b_e_Intercept, b_eh_Intercept, b_th_Intercept) %>% 
    rename(r_tref = b_rtref_Intercept, 
           e = b_e_Intercept,
           eh = b_eh_Intercept,
           th = b_th_Intercept) %>% 
    pivot_longer(everything(), names_to = "parameter") %>% 
    mutate(type = "Posterior")

  prior_draws <- fit_sens %>%
    as_draws_df() %>%
    dplyr::select(prior_b_rtref, prior_b_e, prior_b_eh, prior_b_th) %>% 
    rename(r_tref = prior_b_rtref, 
           e = prior_b_e,
           eh = prior_b_eh,
           th = prior_b_th) %>% 
    pivot_longer(everything(), names_to = "parameter") %>% 
    mutate(type = "Prior")
  
  dist <- bind_rows(prior_draws, post) %>% 
    mutate(iter = i)
  
  dist_list[[i]] <- dist
  
}

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

mutate: new variable 'iter' (integer) with one unique value and 0% NA

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

Warning: There were 13 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

mutate: new variable 'iter' (integer) with one unique value and 0% NA

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

Warning: There were 4 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

mutate: new variable 'iter' (integer) with one unique value and 0% NA

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

mutate: new variable 'iter' (integer) with one unique value and 0% NA

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

mutate: new variable 'iter' (integer) with one unique value and 0% NA

dist_df <- bind_rows(dist_list)

write_csv(dist_df, paste0(home, "/output/prior_sens.csv"))

ggplot(dist_df, aes(value, fill = type)) +
  geom_density(color = NA, alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  ggh4x::facet_grid2(iter~parameter, scales = "free", independent = "y") +
  #coord_cartesian(xlim = c(0, 50)) +
  labs(fill = "") + 
  labs(x = "Value", y = "Density")

ggsave(paste0(home, "/figures/supp/prior_sensi.pdf"), width = 17, height = 17, units = "cm")

ggplot(dist_df %>% filter(parameter == "r_tref"), aes(value, color = as.factor(iter), 
                    linetype = as.factor(type))) +
  geom_density(fill = NA) +
  coord_cartesian(xlim = c(0, 4)) +
  scale_color_viridis(discrete = TRUE) +
  facet_wrap(~parameter, scales = "free") +
  labs(fill = "", linetype = "", color = "iteration")

filter: removed 240,000 rows (75%), 80,000 rows remaining

ggplot(dist_df %>% filter(parameter == "e"), aes(value, color = as.factor(iter), 
                    linetype = as.factor(type))) +
  geom_density(fill = NA) +
  coord_cartesian(xlim = c(0, 12)) +
  scale_color_viridis(discrete = TRUE) +
  facet_wrap(~parameter, scales = "free") +
  labs(fill = "", linetype = "", color = "iteration")

filter: removed 240,000 rows (75%), 80,000 rows remaining

ggplot(dist_df %>% filter(parameter == "eh"), aes(value, color = as.factor(iter), 
                    linetype = as.factor(type))) +
  geom_density(fill = NA) +
  coord_cartesian(xlim = c(0, 10)) +
  scale_color_viridis(discrete = TRUE) +
  facet_wrap(~parameter, scales = "free") +
  labs(fill = "", linetype = "", color = "iteration")

filter: removed 240,000 rows (75%), 80,000 rows remaining

ggplot(dist_df %>% filter(parameter == "th"), aes(value, color = as.factor(iter), 
                    linetype = as.factor(type))) +
  geom_density(fill = NA) +
  coord_cartesian(xlim = c(0, 35)) +
  scale_color_viridis(discrete = TRUE) +
  facet_wrap(~parameter, scales = "free") +
  labs(fill = "", linetype = "", color = "iteration")

filter: removed 240,000 rows (75%), 80,000 rows remaining

dist_df %>% 
  group_by(iter, parameter, type) %>% 
  summarise(mean = mean(value),
            sd = sd(value)) %>% 
  ggplot(aes(iter, mean, color = type)) +
  geom_point(position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0,
                position = position_dodge(width = 0.2)) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free")

group_by: 3 grouping variables (iter, parameter, type)

summarise: now 40 rows and 5 columns, 2 group variables remaining (iter, parameter)

dist_df %>% 
  group_by(iter, parameter, type) %>% 
  summarise(mean = mean(value)) %>% 
  ggplot(aes(iter, mean, color = type)) +
  geom_point(position = position_dodge(width = 0.2)) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free")

group_by: 3 grouping variables (iter, parameter, type)

summarise: now 40 rows and 4 columns, 2 group variables remaining (iter, parameter)

dist_df %>% 
  group_by(iter, parameter, type) %>% 
  summarise(mean = mean(value),
            sd = sd(value)) %>% 
  arrange(parameter) %>% 
  as.data.frame()

group_by: 3 grouping variables (iter, parameter, type)

summarise: now 40 rows and 5 columns, 2 group variables remaining (iter, parameter)

   iter parameter      type       mean         sd
1     1         e Posterior  1.5509848  0.5452377
2     1         e     Prior  1.5473450  1.0585789
3     2         e Posterior  1.8281801  0.6935869
4     2         e     Prior  2.7340863  1.9825645
5     3         e Posterior  1.9445711  0.7507373
6     3         e     Prior  3.9280326  2.8740528
7     4         e Posterior  2.0308057  0.7851491
8     4         e     Prior  5.7526175  4.2679463
9     5         e Posterior  2.0738991  0.8319995
10    5         e     Prior  8.3706061  6.2696777
11    1        eh Posterior  2.2755374  0.3767258
12    1        eh     Prior  2.2667317  1.2801426
13    2        eh Posterior  2.4677415  0.5257818
14    2        eh     Prior  3.3128262  2.1740727
15    3        eh Posterior  2.5565469  0.5719726
16    3        eh     Prior  4.4160862  3.0878564
17    4        eh Posterior  2.6225780  0.6098055
18    4        eh     Prior  6.2134505  4.4571278
19    5        eh Posterior  2.6568042  0.6551016
20    5        eh     Prior  8.8710538  6.5123861
21    1    r_tref Posterior  0.5584244  0.2401727
22    1    r_tref     Prior  0.7241750  0.5076331
23    2    r_tref Posterior  0.7053907  0.3803210
24    2    r_tref     Prior  1.0239639  0.7458058
25    3    r_tref Posterior  0.7840347  0.4564422
26    3    r_tref     Prior  1.4564463  1.0584843
27    4    r_tref Posterior  0.8389461  0.5159858
28    4    r_tref     Prior  2.1424920  1.5787953
29    5    r_tref Posterior  0.8858437  0.5950614
30    5    r_tref     Prior  3.1128160  2.3175503
31    1        th Posterior  7.3090808  2.1976232
32    1        th     Prior 11.2071689  6.3963478
33    2        th Posterior  6.5898775  2.1027913
34    2        th     Prior 13.7260271  8.6061495
35    3        th Posterior  6.3742811  2.1228914
36    3        th     Prior 17.9937275 12.0153405
37    4        th Posterior  6.2053125  2.0519829
38    4        th     Prior 24.2390235 16.9845314
39    5        th Posterior  6.1379047  2.0948903
40    5        th     Prior 34.5026161 24.9471179

Calculate T_opt

ts <- dat %>% 
  data_grid(temp = seq_range(temp, n = 100)) %>% 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) %>% 
  add_epred_draws(fit_global, re_formula = NA, ndraws = 8000) 

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

p1 <- ggplot(ts) + 
  geom_line(aes(temp, y = .epred, group = .draw), alpha = .1) + 
  coord_cartesian(ylim = c(0.1, 0.4))

p2 <- ggplot(ts) + 
  geom_line(aes(temp, y = .epred, group = .draw), alpha = .1) +
  coord_cartesian(ylim = c(0.1, 0.4))

p1 + p2

# Compute quantiles across spaghetties
t_opt <- ts %>% 
  group_by(.draw) %>% 
  filter(.epred == max(.epred)) %>% 
  ungroup() %>% 
  filter(!temp == min(temp)) %>% # remove values where there was no optimum
  filter(!temp == max(temp))

group_by: one grouping variable (.draw)

filter (grouped): removed 792,000 rows (99%), 8,000 rows remaining

ungroup: no grouping variables

filter: removed 2 rows (<1%), 7,998 rows remaining

filter: removed one row (<1%), 7,997 rows remaining

# These are used in the fit-temp script
quantile(t_opt$temp, probs = c(0.025, 0.5, 0.975))

     2.5%       50%     97.5% 
 8.771866  9.468607 10.304697 

ggplot(t_opt, aes(temp)) + 
  geom_histogram(bins = length(unique(t_opt$temp))) + 
  theme_void()

# Make the main plot (conditional effect of temperature, with and without random effects)
# Predictions without random effects
glob <- dat %>%
  data_grid(temp = seq_range(temp, n = 100)) %>% 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) %>% 
  add_epred_draws(fit_global, re_formula = NA) %>% #
  ungroup()

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

dat %>%
  group_by(area) %>%
  data_grid(temp = seq_range(temp, n = 100)) %>% 
  ungroup() %>% 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) %>% 
  add_epred_draws(fitbs) %>% 
  ungroup() %>% 
  ggplot(aes(temp, y = .epred, color = factor(area, order$area))) +
  stat_lineribbon(data = glob, aes(temp, .epred), .width = c(0.9), inherit.aes = FALSE,
                 fill = "black", color = "black", alpha = 0.1) +
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 1, alpha = 0.8) +
  stat_lineribbon(.width = c(0)) +
  stat_lineribbon(data = glob, aes(temp, .epred), .width = c(0), inherit.aes = FALSE,
                 color = "black", alpha = 0.9, linetype = 2) +
  coord_cartesian(expand = 0) +
  scale_color_manual(values = colors, name = "Site") +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)") + 
  NULL

group_by: one grouping variable (area)
ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

ggsave(paste0(home, "/figures/sharpe_school_bayes.pdf"), width = 17, height = 17, units = "cm")
#ggsave(paste0(home, "/figures/for-talks/sharpe_school_bayes.pdf"), width = 14, height = 14, units = "cm")

Area-specific T_opts as a ridgeplot

# Calculate T_opt by area
tsa <- dat %>% 
  data_grid(area = unique(dat$area),
            temp = seq_range(temp, n = 100)) %>% 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) %>% 
  add_epred_draws(fitbs, re_formula = NULL, ndraws = 10000) 

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

t_opt_a <- tsa %>% 
  group_by(.draw, area) %>% 
  filter(.epred == max(.epred)) %>% 
  ungroup() %>% 
  filter(!temp == min(temp)) %>% # remove values where there was no optimum
  filter(!temp == max(temp))

group_by: 2 grouping variables (.draw, area)

filter (grouped): removed 9,900,000 rows (99%), 100,000 rows remaining

ungroup: no grouping variables

filter: removed 4,179 rows (4%), 95,821 rows remaining

filter: removed 8,461 rows (9%), 87,360 rows remaining

t_opt_a_sum <- t_opt_a %>% 
  group_by(area) %>% 
  summarise(median = median(temp),
            lwr = quantile(temp, probs = 0.1),
            upr = quantile(temp, probs = 0.9))

group_by: one grouping variable (area)

summarise: now 10 rows and 4 columns, ungrouped

# How many datapoints are below or above optimum by area and globally?
left_join(dat, t_opt_a_sum, by = "area") %>% 
  filter(temp < median) %>% 
  nrow()

left_join: added 3 columns (median, lwr, upr)

           > rows only in x     0

           > rows only in y  (  0)

           > matched rows     306

           >                 =====

           > rows total       306

filter: removed 117 rows (38%), 189 rows remaining

[1] 189

left_join(dat, t_opt_a_sum, by = "area") %>% 
  mutate(above = ifelse(temp > median, "Y", "N")) %>% 
  group_by(area) %>% 
  count(above) %>% 
  pivot_wider(names_from = above, values_from = n) %>% 
  mutate(prop = Y/(Y+N)) %>% 
  arrange(desc(prop))

left_join: added 3 columns (median, lwr, upr)

           > rows only in x     0

           > rows only in y  (  0)

           > matched rows     306

           >                 =====

           > rows total       306

mutate: new variable 'above' (character) with 2 unique values and 0% NA

group_by: one grouping variable (area)

count: now 16 rows and 3 columns, one group variable remaining (area)

pivot_wider: reorganized (above, n) into (N, Y) [was 16x3, now 10x3]

mutate (grouped): new variable 'prop' (double) with 7 unique values and 40% NA

# A tibble: 10 × 4
# Groups:   area [10]
   area      N     Y    prop
   <chr> <int> <int>   <dbl>
 1 FM       10    27  0.730 
 2 BT        7    15  0.682 
 3 SI_EK    14    10  0.417 
 4 JM       37    23  0.383 
 5 FB       44     3  0.0638
 6 BS       16     1  0.0588
 7 HO       29    NA NA     
 8 MU       18    NA NA     
 9 RA       14    NA NA     
10 SI_HA    NA    38 NA     

# Without impacted?
left_join(dat, t_opt_a_sum, by = "area") %>% 
  filter(!area %in% c("BT", "SI_EK")) %>% 
  filter(temp < median) %>% 
  nrow()

left_join: added 3 columns (median, lwr, upr)

           > rows only in x     0

           > rows only in y  (  0)

           > matched rows     306

           >                 =====

           > rows total       306

filter: removed 46 rows (15%), 260 rows remaining

filter: removed 92 rows (35%), 168 rows remaining

[1] 168

dat %>% 
  group_by(area) %>% 
  summarise(n = n()) %>% 
  arrange(desc(n))

group_by: one grouping variable (area)

summarise: now 10 rows and 2 columns, ungrouped

# A tibble: 10 × 2
   area      n
   <chr> <int>
 1 JM       60
 2 FB       47
 3 SI_HA    38
 4 FM       37
 5 HO       29
 6 SI_EK    24
 7 BT       22
 8 MU       18
 9 BS       17
10 RA       14

# nrow(filter(dat, temp < 9.468607))
# nrow(dat)

ggplot(t_opt, aes(temp)) + 
  geom_histogram(bins = length(unique(t_opt$temp)))

# Compute quantiles across spaghetties
rect_dat <- data.frame(xmin = c(quantile(t_opt$temp, probs = c(0.05)), quantile(t_opt$temp, probs = c(0.1))),
                       xmax = c(quantile(t_opt$temp, probs = c(0.95)), quantile(t_opt$temp, probs = c(0.9))),
                       ymin = c(-Inf, -Inf),
                       ymax = c(Inf, Inf),
                       interval = c("90%", "80%"))

ridge <- ggplot(t_opt_a, aes(temp, factor(area, levels = rev(order$area)), fill = factor(area, levels = rev(order$area)))) + 
  geom_vline(xintercept = quantile(t_opt$temp, probs = c(0.5)), alpha = 0.6, linetype = 2, linewidth = 0.75) +
  geom_rect(data = rect_dat %>% filter(interval == "90%"), aes(xmin = xmin, xmax = xmax), ymin = -Inf, ymax = Inf, 
            inherit.aes = FALSE, fill = "grey", alpha = 0.35) +
  geom_rect(data = rect_dat %>% filter(interval == "80%"), aes(xmin = xmin, xmax = xmax), ymin = -Inf, ymax = Inf, 
            inherit.aes = FALSE, fill = "grey", alpha = 1) +
  geom_density_ridges(alpha = 0.7, color = NA, scale = 1,
                      stat = "binline", bins = 50) +
  scale_fill_manual(values = rev(colors), name = "Site") +
  scale_color_manual(values = rev(colors), name = "Site") +
  labs(x = "Temperature (°C)", y = "") +
  guides(fill = "none") +
  geom_errorbar(data = t_opt_a_sum, aes(xmin = lwr, xmax = upr, y = area, color = factor(area, levels = rev(order$area))), linewidth = 1,
                inherit.aes = FALSE, width = 0) +
  geom_point(data = t_opt_a_sum, aes(median, area, fill = factor(area, levels = rev(order$area))), shape = 21, color = "white", size = 3) +
  NULL

filter: removed one row (50%), one row remaining

filter: removed one row (50%), one row remaining

ggsave(paste0(home, "/figures/t_opt_ridges.pdf" ), width = 17, height = 17, units = "cm")

# Any patterns in the residuals (area spec and global)
t_opt_a_sum$global <- quantile(t_opt$temp, probs = c(0.5))

t_opt_a_sum$diff <- t_opt_a_sum$median - t_opt_a_sum$global

order2 <- VBG %>%
  ungroup() %>%
  group_by(area) %>%
  summarise(mean_temp = mean(temp)) %>%
  arrange(desc(mean_temp))

ungroup: no grouping variables

group_by: one grouping variable (area)

summarise: now 10 rows and 2 columns, ungrouped

t_opt_a_sum <- left_join(t_opt_a_sum, order2)

Joining with `by = join_by(area)`

left_join: added one column (mean_temp)

           > rows only in x    0

           > rows only in y  ( 0)

           > matched rows     10

           >                 ====

           > rows total       10

diff <- ggplot(t_opt_a_sum, aes(mean_temp, diff, label = area)) +
  geom_point(color = "gray10") + 
  geom_errorbar(aes(ymin = lwr-global, ymax = upr-global), width = 0, alpha = 0.4) +
  labs(x = "Mean site temperature (°C)", y = "Site-specific - Global optimum") +
  geom_hline(yintercept = 0, alpha = 0.5, color = "gray30", linetype = 2) +
  NULL
  
(ridge & coord_flip() & guides(color = "none")) + diff + plot_layout(widths = c(1, 0.6)) + plot_annotation(tag_levels = "A")

ggsave(paste0(home, "/figures/t_opt_ridges.pdf" ), width = 17, height = 8, units = "cm", device = cairo_pdf)

Plot area specific predictions

# Compare with area-specific sharpe scool!
area_pred_brms <- dat %>%
  group_by(area) %>% 
  data_grid(temp = seq_range(temp, n = 100)) %>% 
  ungroup() %>% 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) %>% 
  add_epred_draws(fitbs)

group_by: one grouping variable (area)

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

# Plot area specific predictions as facets
p1 <- ggplot(area_pred_brms, aes(temp, y = .epred, color = factor(area, order$area),
                           fill = factor(area, order$area))) +
  stat_lineribbon(.width = c(0.95), alpha = 0.4, color = NA) +
  stat_lineribbon(.width = c(0)) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Site") +
  scale_fill_manual(values = colors, name = "Site") +
  scale_linetype_manual(values = 2) +
  facet_wrap(~factor(area, rev(order$area))) +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)",
       linetype = "")

p1

ggsave(paste0(home, "/figures/supp/sharpe_school_bayes_ci_facet.pdf" ), width = 17, height = 17, units = "cm")
ggsave(paste0(home, "/figures/for-talks/sharpe_school_bayes_ci_facet.pdf" ), width = 14, height = 14, units = "cm")